######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure G1
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#########################
#Estimation for Figure G.1
#########################
#----------------
#(1) Load result
#----------------
load("result/bpcausal_original.RData")

#Summary
sout1 = coefSummary(result.original)

#--------------
#(2) Parameters
#--------------
#Parameters
newx = 1992:2019
cex = 4
varname = c("GDP per capita (log)",
            "Population (log)",
            "Election",
            "FDI inflow (% GDP)",
            "Debt service (% GNI)",
            "ODA received (% GNI)",
            "Polity",
            "Temporary UNSC member",
            "UNGA Voting (Ideal Point distance from US)")

#########################
#Produce Figure G1
#########################
#Setup
png(filename = "figure/Figure_G1.png",
    height = 2400,
    width = 3000)

par(mfrow = c(3, 3),
    mar = c(4, 7, 5, 1),
    lend = 1)

#Plot
for (i in 1:length(varname)){
  #data
  out = sout1$est.xi[[i+1]]
  ymin = floor(min(out$ci_l)*10) / 10
  ymax = ceiling(max(out$ci_u)*10) / 10
  step = (ymax - ymin)/4
  step = ifelse(i == 5, 0.2, step)
  
  ylim = c(ymin, ymax)
  
  #empty plot
  plot(1, 
       type = "n",
       xlab = "",
       ylab = "",
       axes = F,
       xlim = range(newx),
       ylim = ylim)
  box()
  #Title
  title(main = varname[i],
        line = 1,
        cex.main = cex)
  #X-axis
  axis(side = 1,
       at = seq(min(newx),
                max(newx),
                by = 2),
       cex.axis = cex)
  #Y-axis
  axis(side = 2,
       at = seq(min(ylim),
                max(ylim),
                by = step),
       cex.axis = cex)
  mtext(text = "Coefficients",
        side = 2,
        line = 3,
        cex = cex-1)
  #Ablines
  abline(h = 0,
         col = "gray20",
         lty = 2,
         lwd = cex)
  
  #Mean
  lines(x = newx,
        y = out$mean,
        col = 1,
        lty = 1,
        lwd = cex)
  polygon(x = c(rev(newx),
                newx),
          y = c(rev(out$ci_l),
                out$ci_u),
          col = "#55555530",
          border = NA)
}
#Close
dev.off()